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Numerical lattice gauge theory computations to generate gauge field configurations including the 
effects of dynamical fermions are usually carried out using algorithms that require the molecu- 
lar dynamics evolution of gauge fields using symplectic integrators. Sophisticated integrators are 
in common use but are hard to optimise, and force-gradient integrators show promise especially 
for large lattice volumes. We explain why symplectic integrators lead to very efficient Monte 
Carlo algorithms because they exactly conserve a shadow Hamiltonian. The shadow Hamiltonian 
may be expanded in terms of Poisson brackets, and can be used to optimize the integrators. We 
show how this may be done for gauge theories by extending the formulation of Hamiltonian me- 
chanics on Lie groups to include Poisson brackets and shadows, and by giving a general method 
for the practical computation of forces, force-gradients, and Poisson brackets for gauge theories. 
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I. INTRODUCTION 

Essentially all algorithms used in lattice gauge theory 
computations to generate gauge field configurations in- 
cluding the effects of dynamical fermions are variants 
of the Hybrid Monte Carlo (HMC) algorithm 0, which 
requires a reversible and area-preserving integrator for 
its molecular dynamics step. The simplest such method 
is the leapfrog integrator, but there is a large class of 
symplectic integrators that have these properties and 
are potentially more cost-effective. Indeed, many state- 
of-the-art computations use the second order minimum 
norm integrator [^-^ which has a free parameter, which 
heretofore has been tuned in an ad hoc manner. 

The formulation of Hamiltonian dynamics on Lie group 
manifolds, which is required for molecular dynamics on 
gauge fields ^, and the fact that symplectic integra- 
tors conserve a shadow Hamiltonian are well known; our 
goal is to combine the two and show how to construct the 
shadow Hamiltonian for gauge theories. This is most eas- 
ily done using the formalism of differential forms [p|-p2[; 
in order to fix our notation and establish the necessary 
results, some of which are not easy to find in the litera- 
ture, we provide a brief review in Appendix 

The shadow Hamiltonian is expressed as an asymp- 
totic expansion in the integration step size St whose co- 
efficients depend on the parameters specifying the inte- 
grator under consideration and a collection of Poisson 
brackets. These Poisson brackets are complicated func- 



tions on phase space, where in the case of gauge field 
molecular dynamics a point in phase space is an en- 
tire gauge field configuration and its associated "ficti- 
tious" momenta. For extensive systems such as field the- 
ories, unlike the few body systems considered previously 
Q , the values of the Poisson brackets have a distri- 
bution that is sharply peaked about their mean values 
when we choose the starting points of their molecular 
dynamics trajectories to be chosen from the distribution 
e~^ , as is done in the HMC algorithm. This may be un- 
derstood as a consequence of the central limit theorem 
applied to the contributions to the Poisson brackets com- 
ing from many independent regions of space-time. This 
means that for configurations that occur with non van- 
ishingly small probability the shadow Hamiltonian may 
be considered to be a function of the average values of 
the Poisson brackets; if these are measured on a few test 
trajectories then the integrator parameters may be cho- 
sen to minimize the computational cost |15-17. Per- 



haps surprisingly this does not correspond to minimizing 
the average difference between the Hamiltonian and its 
shadow 1^ and instead to minimizing the variance of the 
distribution of the shadow. We shall not be concerned 
with the details of this tuning procedure here, but we 
refer the interested reader to Ijl^l for details: instead, the 



^ Since the shadow is only defined up to an additive constant this 
cannot be too surprising. 
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aim of this paper is to explain how the Poisson brack- 
ets, forces, and force-gradients may be computed at any 
given point in phase space. 

In |j expressions for the molecular dynamics force 
were derived from the classical mechanics specified by 
the Hamiltonian function and a suitable chosen group- 
invariant fundamential two-form. We extend this anal- 
ysis to obtain an expression for the force-gradient for 
gauge fields [18 , which can be used to provide a "sec- 
ond derivative" integrator step for the construction of 



improved integrators |13, 14 



A. Multiple link updates 

For much of this paper we shall be considering a Hamil- 
tonian system with a phase space which is the cotangent 
bundle T*Q over a base space that is a Lie group mani- 
fold Q and whose fibres are isomorphic to its Lie algebra. 
We shall call the cotangent one-forms "momenta", al- 
though in the context of HMC they are called "fictitious 
momenta" as they are quite different from the canoni- 
cal momenta of the underlying field theory. For a gauge 
field theory we may associate such a phase space with 
every link of the lattice. One might at first think that 
we need to introduce some fibre bundle structure over the 
space-time lattice itself, but fortunately that is not neces- 
sary. We can consider the molecular dynamics evolution 
of each gauge link separately; they are coupled together 
through the action that plays the role of the potential 
energy part of the Hamiltonian, but the kinetic energy 
part does not couple different links. For HMC we are 
free to choose the form of the kinetic energy, so we can 
take it to be of the form ^ T{p) = i J2i (^iPi where pi 
is the momentum associated with the link i, and cg is a 
link-dependent coefficient that is constant in molecular 
dynamics "fictitious" time. If we wish to evolve the sin- 
gle link (f)i on its own we can choose ce> — bi±i so that 
= dUjdvf = dT/dpt = crpr = if ^ 7^ We 
are also free to choose cg — 1 for all links, which is the 
usual situation where we update the gauge field simul- 
taneously across the entire lattice. Another interesting 
choice for the kinetic energy is to choose ci = 1 for all 
spatial links and Q = ^ for all temporal ones: this is the 
procedure suggested in [|l9[ |^ for evolving anisotropic 
lattices ^. The momentum anisotropy ^ is a parameter 
that can be adjusted to optimize the HMC algorithm 
for a given anisotropy in the action; if the spatial and 
temporal contributions to the Poisson brackets are mea- 
sured separately then the techniques of 17 can be used 



to tune ^ along with other integrator parameters. 



^ For notational simplicity we consider here a theory with a scalar 
field and the corresponding momentum p defined on the links 
of a lattice. 

^ In po[ the temporal step size is adjusted rather than the kinetic 
energy, but this is equivalent after a rescaling of the temporal 
momenta. 



B. Pseudofermion forces 

So far we have only been discussing pure gauge theo- 
ries, but in practice the cost of most lattice computations 
is dominated by the inclusion of fermions. This is because 
we need to solve a large system of linear equations in or- 
der to update the fictitious moments (i.e., to apply the 
Hamiltonian vector field S in the notation we will intro- 
duce later). Typically we have an action S which is the 
sum of a pure gauge part Sg, built out of sums of small 
Wilson loops (traces of a closed loops of gauge links) such 
as plaquettes, and a pseudofermion part Sp built out of 
sums of pairs of pseudofermion fields (j) connected by a 
string of gauge links. If we want to compute the force 
acting on a particular gauge link ^ U then it is conve- 
nient to write Sg = Retr(nC/) and Sp = (f)^ M~'^ {U)(t) 
where the "staple" □ is the sum of all gauge link strings 
that connect the ends of the link U that correspond to 
the Wilson loops in Sg , and the Hermitian lattice matrix 
M.{U) is the sum of all gauge link strings that include U 
that occur in Sp- For a local action all of these strings 
are in some neighbourhood of J7, and we have dropped all 
other terms in the action because they are independent 
of U and therefore do not contribute to the force on that 
particular link. In reality we update many or all the links 
on the lattice at once, so we compute the force on each 
link in parallel. By the "force" we mean the quantity 
ej(S')T* where is a linear differential operator (vector 
field) whose action on U is specified by e,:(i7) = —TiU 
and which we shall define carefully later (|l^), and is 
the representation of a generator of the gauge group. It 
is important to note that here Bi acts only on the gauge 
link [/, it gives zero if applied to any other link variable. 
There is an opportunity for confusion when we refer to 
Bi as a vector field; it is a vector field defined over the 
phase space of the link U , but it is not a field over the 
space-time lattice. In order to reduce confusion we refer 
to quantities defined over the space-time lattice as lattice 
vectors^ and space time linear differential operators such 
as the Dirac operator (or more precisely lattice difference 
operators acting on lattice vectors such as the Wilson- 
Dirac operator) as lattice matrices. 

The contribution to the force from the pure gauge 
part of the action is ei{SG)T' = Retr(^nei([/))r' = 

-Retr(nT,C/)T* = - Retr(f7nr,)T' = -ar(C/n), 
T being the projector onto the Lie algebra, that is 
T{X) = Retr(XTi)rYo where there is an implicit sum 
over i as usual and the generators Ti are normalized such 
that tr{TiTj) — aSij. If the gauge group is SU(iV) and 
we choose its generators to be anti-Hermitian so as not 
to introduce artificial factors of i, then T{X) is just the 
traceless anti-Hermitian part of X. 



* We shall refer both to a gauge link variable and the link on which 
it lives as U when there is no ambiguity. 
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The pseudofermion contribution to the 
e,{SF)T' = 0te,(x([/)-i)0r\ Since 



ear differential 



= -M 

'■r(M{U) 



operator we have 

e,{M)M-^ + Me,{M 
-1 



-1 



force is 

Bi is a hn- 

= e,(I) = 
and hence 
Therefore ei{SF)T^ = 

where we have defined 



tion of the group. Moreover, the fact that Hamiltonian 
vector fields form a Lie algebra is crucial for the defini- 



-Retr e,(^M{U)jX d) 

X = M.~^ (j) to he the solution of a large but sparse sys- 
tem of linear equations (since M is local), this may be 
computed on all lattice sites and used to update some or 
all gauge links in parallel. The outer product X ® X'^ is 
the rank one Hermitian lattice matrix whose action on an 
arbitrary lattice vector y is proportional to the projection 
of y along X, namely [X ® X'^)y ^ X{X'^y). 

We can express the pseudofermion action in the form 
Sf = - Retr[A4(t/)Xi8Xt] analogous to that of Sg if we 
consider X to be a lattice vector that is independent of U. 
This means that once we have computed X the calcula- 
tion of the gauge and pseudofermion parts of the force 
and related quantities are very similar. Both the gauge 
and pseudofermion actions can be written as the trace of 
lattice operators times U, where the lattice operators are 
either local (n and A4) or low rank {X®X'^). Both local 
and low rank operators are relatively cheap to apply to 
lattice vectors or to trace, the former only involving links 
in the neighbourhood of U, and the latter only involving 
inner products of lattice vectors. For example, we may 
evaluate the trace tr[M{U)X iS> X'f] = X''M{U)X as the 
inner product of X'' with the vector j\4{U)X. 

If we include spin degrees of freedom then we must 
replace X (3 X'^ by a sum of outer products for each 
spin component, but the result is still a low rank ma- 
trix which is therefore cheap to apply. Likewise if we 
wish to introduce n pseudofermion fields so as to reduce 
the noise in the stochastic estimate of the fermionic force 
and thus defer the breakdown in the asymptotic expan- 
sion for the shadow Hamiltonian to significantly larger 
integrator step sizes [pTj-psI, then we only increase the 
rank by a factor of n. 



C. Outline 

The structure of this paper is as follows. In §0 we con- 
sider the general formulation of Hamiltonian mechanics 
on a symplectic manifold this serves to introduce 
the important concepts of the fundamental 2-form, the 
Hamiltonian vector field it associates with any 0-form, 
and the Poisson bracket of two 0-forms. We show that 
Poisson brackets satisfy the Jacobi identity, and that the 
commutator of two Hamiltonian vector fields is itself a 
Hamiltonian vector field, and explain the isomorphism 
between the Lie algebra of commutators of Hamiltonian 
vector fields and that of Poisson brackets of 0-forms. The 
reason we need all this mathematical machinery is that 
when we consider Hamiltonian mechanics on Lie groups 
in § IV we will introduce a non-trivial fundamental 2-form 



tion of the shadow Hamiltonian, which we give in §IH 



The exposition assumes some knowledge of the theory of 
differential forms, an overview of which is given in Ap- 
pendix 



§III| introduces symplectic integrators by noting that if 
a 0-form on phase space only depends on the momenta p 
or only on the positions q then the integral curves of its 
Hamiltonian vector field are easily found. We are inter- 
ested in Hamiltonians H{q,p) — T{p) + S{q) that are the 
sum of two such functions, and we show how this allows 
us to construct symplectic integrators to find approxi- 
mate integral curves for H using the Baker-Campbell- 
Hausdorff (BCH) formula. We give some simple examples 
of integrators for a system on a symplectic manifold with 
fundamental 2-form oj — dq /\ dp, and show how to com- 
pute the corresponding shadow Hamiltonians. When the 
kinetic energy is of the form Tijj) — ip^ we show that 
the Poisson bracket {5*, {S,T}} is independent of p and 
explain how it may thus be used to construct a force- 
gradient integrator step. 



f:IV defines a symplectic structure on Lie group man- 
ifolds, or more precisely on their cotangent bundle T*Q, 
that is compatible with the group structure. This is done 
by introducing the natural fundamental 2-form terms of 
Maurer-Cartan forms, and it is here that the mathe- 
matical framework we have developed becomes necessary. 
We derive explicit formula for Hamiltonian vector fields 
and Poisson brackets in terms of the momentum coor- 
dinates (which are well-defined globally) and the family 
of left-invariant vector fields dual to the Maurer-Cartan 
forms. All the independent Poisson brackets of S and 
T that can occur in shadow Hamiltonians up to and in- 
cluding 0{5t'^) are given explicitly for the case where S 
is momentum-independent and T is quadratic in the mo- 
menta. We then show how to express the results in terms 
of matrix representations of the Lie group, as these are 
what is used in practice. 

In §^ we evaluate the formulae for the Poisson brack- 
ets for the physically interesting case of the fundamental 
representation of SU(A^). We show that they can all be 
expressed as traces of a collection of Lie-algebra- valued 
quantities: as these live on links we name them basic 
lattice vectors. 



in order to make the dynamics symmetric under the ac- 



In § VI we address the problem of computing these ba- 
sic lattice vectors. We do this first for the simple case 
where only a single link is updated, and then introduce 
the algebra of towers to give an efficient way of comput- 
ing them in general. 

Appendix ^ gives a brief survey of the theory of differ- 
ential forms and serves to fix our notation and conven- 
tions, as does Appendix |^ which gives an overview of the 
properties of Lie groups. 
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II. HAMILTONIAN MECHANICS 



and also 



A. Symplectic Manifolds 

A Hamiltonian system is defined on phase space which 
is a differential manifold Ai with a symplectic struc- 
ture given by some fundamental 2-form uj that is closed, 
du} — 0, and globally invertible. Phase space is usually 
the cotangent bundle T*Q over some configuration space 
manifold Q. For every 0-form F G A° on A^, that is for 
every C°° smooth function F : M. ^ there is a corre- 
sponding Hamiltonian vector field F G Ham M. such that 
dF = ipijL): in other words dF{y) — {ipU)){y) = lj{F, y) 
for any vector field y. 

A 0-form Z corresponds to a vanishing Hamiltonian 
vector field iff dZ — 0, so we have the following short 
exact sequence ^> IR AP{M) -> HamA^ — s> 0. 
This implies that there is a bijective diffeomorphism 
A°(7W)/R o HamTU. The nature of this correspon- 
dence between 0-fornis (up to an additive constant) and 
Hamiltonian vector fields will be examined further in the 
following sections. 



B. Poisson Brackets 

Consider the action of a Hamiltonian vector field F on 
a 0-form G, 

FG = dG{F) = i^u{F) = lj{G, F) = {F, G}, 

where in the first equality we have made use of the def- 
inition of the exterior derivative of a 0-form G acting 
on an arbitrary vector field y, dG{y) = yG, and in 
the last equality we have introduced the Poisson bracket 
{A,B} = —iji}{A,B) for any pair of 0-forms A and B. 
The minus sign has to appear somewhere, and our con- 
vention is to introduce it here in the definition of the 
Poisson bracket. 



C. Jacobi Identity 



The invariant expression (A3) for the exterior deriva- 
tive doj of a 2-form ui applied to three arbitrary vector 
fields X, y, and z 

duj{x, y, z) = xu}{y, z) + yu:{z, x) + zu}{x, y) 

-Lo{[x,y\,z) - u}{[y,z\,x) - u:{[z,x\,y), 

displays an interesting cyclic symmetry in the three vec- 
tor fields x,y, and z. This has an important consequence 
if a; is the fundamental 2-form and the vector fields are 
Hamiltonian; if A^ B, and G are three arbitrary 0-forms 
then 

Alj{B, C) = -A{B, G} = -{A, {B, C}}, 



a;([l, B], C) = -u{C, [A, B]) = -dC([A, B]) 
= -[A, B]G = [BA - AB)G 

We thus find that the condition du) — implies that 
the cyclic sum of of nested Poisson brackets must 
vanish, du{A,B,C) = {A,{B,G}} + {B,{G,A}} + 
{C, {A, _B}} = 0: this is just the Jacobi identity which, 
together with the antisymmetry of the Poisson bracket, 
demonstrates that 0-forms on Ai together with the prod- 
uct given by the Poisson bracket form a Lie algebra. 

We can use the Jacobi identity to derive another useful 
result. The commutator of any two vector fields is a 
vector field (q.v., equation (A2)); if both vector fields are 



Hamiltonian then their commutator is also Hamiltonian, 
since 

[A, B]G = {AB - BA)G = {A, {B, G}} - {B, {A, G}} 
= -{G, {A, B}} = {{A, B}, G} = {ASyG 

where we applied the Jacobi identity in the antepenulti- 
mate step. Since this must hold VC £ A*' we have 



[A,B] = {A,B} G HamA^ 



(1) 



telling us that not only is the commutator of two Hamil- 
tonian vector fields Hamiltonian as promised, but also 
that it corresponds to the 0-form that is the Poisson 
bracket of the 0-forms corresponding to the original pair 
of Hamiltonian vector fields. The bijection A°(7V()/R o 
Ham(A^) is therefore an isomorphism of Lie algebras. 



D. Lie Derivatives and Equations of Motion 

Given a Hamiltonian H € A°(A^) and a fundamental 
2-forni Lij we may construct the Hamiltonian vector field 
H, and for any point p ^ A4 we may — at least locally 
— define an integral curve. We may also define a local 
flow a :X xU ^ M. oi trajectories starting at any point 
p e C in some neighbourhood of p, a : IR ^ A1, 
satisfying Hamilton's equations da/dt = H and the initial 
condition cr(0) — p. Hamilton's equations are thus most 
naturally expressed in terms of Lie derivatives ( ^A 
dT/dt — CjjT, for any tensor T. In particular a scalar 
field (0-form) F, vector field v, and 1-form 6 must obey 

dF 

— =C^F = HF = {H,F}, 
do 

and — = CfjO = {ifjd + dijj)6. 

The formal solution of the equation of motion 
dT/dt = CfjT is T{t) = exp(t£^)r(0), where 
the exponential function is defined as exp(t£^) ~ 

iim„^^ (1 + ^/:^)" = Y.T=o{tc^yi^'- 
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III. SYMPLECTIC INTEGRATORS AND 
SHADOW HAMILTONIANS 

A. Baker— Campbell— HausdorfF Formula 

The BCH formula states that if A and B belong to an 
associative algebra then 



their exponential maps t M> [exp{tA/n) exp{tB /n)]'^ for 
some n € N. Such a map is called a symplectic integrator 
as it manifestly preserves the symplectic structure since 
each individual exponential map does. The BCH formula 
(||) tells us that this curve is in fact itself the integral 
curve of a vector field Df /„ 



(2) 



where the c„, belonging to the /ree Lie algebra are 
recursively determined from the relations ci = A + B 
and 



{n + l)c„+i 

L"/2J 



Bo 



'-^ (2m)! 

1=1 ^ ' fci 



ad Cfcj ... ad c/j^^ (A + i?) 



fc2m >1 
''l + --- + ''2m = i 



-i(adc„)(^-B) forn>l, (3) 

where ad a : 6 H> [a, 6] and the Bernouilli numbers Bn 
are defined by 



The first few terms in the Hausdorff series are 

ln(e^e-S) = (A+B) + i[A,B] 

+ ^{IA,[A,B]]-[B,[A.B]])-^IB.IA,[A.B]]] 

-A[B^[A,IA,IA,B]]]] -6llA,B],lA,lA,B]]] 
+ 7311 I +4.1B,IB,[A,IA,B]]]] -2llA.B],lB.,lA,B]]] \ + 
-IA,IA,IA,IA.B]]]] +[B,[B,[B,[A,B]]]] 

From this we easily obtain the corresponding formula for 
a symmetric product 

32[B,[B,[A,[A,S]]]] -ie[[A,B],[B,[A.B]]] ^ 
+28[B,[A,[A,[A,B]]]] +12[[A,B],[A,[A,B]]] | +■■■ 
+8[B,[B,[B,[A,B]]]] +7[A,[A,[A,[A,B]]]] 



B. Symplectic Integrators 

The integral curve of a Hamiltonian vector field A is 
given by the exponential map t i— exjp{tA) acting on 
the initial point. Given two Hamiltonian vector fields 
A and B we can construct a curve that is alternately 
tangential to each vector field from the composition of 



^ That is the Lie algebra whose Lie bracket is the commutator 
constructed from the associative product. For more details about 
free Lie algebras and a proof of the BCH formula see Appendix B 
of [0. 



tA 



tB 



exp — exp — 
\ n I \ n 



exp + + Vc,„(A,S)(- 

\ m=2 ^ 

(^A + B+f2^c„,{A,B)^^y 



exp 
exp ( A/n^ 



where = A + B + YZ=2 Cm(A, B)e^''-\ As aU the 
Cm are commutators, equation (|l|) tells us that is a 
Hamiltonian vector field corresponding to the shadow 0- 
form Dg under the isomorphism HamA^ o Ao(A^)/R 
discussed before. In other words, — where the 
0-form D, = A + B + Y.^^2 Cm(^, B)e"'-^ with the c'^ 
defined by (||) in terms of the Poisson bracket image of the 
adjoint under the Lie algebra isomorphism (|^) a.dA 
ad A where ad A : B i-^ {A, B}. We note in passing that 
the shadow is only defined up to an additive constant. 

The BCH formula is obtained by formal manipulation 
of the exponential series, so we should choose a suffi- 
ciently large n to ensure that the Hausdorff series con- 
verges. In order to study the convergence of the BCH for- 
mula we need to specify a topology on the space of Hamil- 
tonian vector fields Ham A4. It is simpler to ask the same 
question about the convergence of the corresponding ex- 
pansion for the shadow Hamiltonian, for which there is 
an obvious topology as the coefficients are 0-forms and 
we can use the usual Lp norms. In most cases of interest 
none of these norms are bounded, so the series is only 
asymptotic at best. In HMC the momenta are selected 
from a Gaussian distribution e~'^^P\ so the values of the 
Poisson brackets can become arbitrarily large, but with 
exponentially small probability. There is no value of s for 
which the Hausdorff series always converges, but it might 
well be that for any S > we can find an e > such that 
it does converge with probability > 1 — S. This may 
be acceptable for HMC, where an exponentially small 
chance of a trajectory becoming unstable is unimportant: 
it will presumably be rejected and the next momentum 
or pseudofermion refreshment will resolve the problem. If 
the large norm comes from the gauge field configuration 
then there could be more severe problems. 
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C. Symmetric Symplectic Integrators 



or equivalently 



In general a symplectic integrator is not reversible, that 
is the group commutator 

exp{--tA/n) exp{—tB /n) cxp{tA/n) exp{tB /n) ^ I; 

indeed we immediately see from this expression that that 
the integrator is reversible iff [A, B] = 0. This blemish 
is easily eradicated by using a symmetric symplectic in- 
tegrator, such as exp{tA/2n) exp{tB /n) exp{tA/2n). An 
additional advantage of such integrators is that only even 
powers of e occur in the Hausdorff series for their shadow 
Hamiltonians D^, so A + B — — 0{e'^), making them 
better approximations to the exponential map of A + B 
itself. 



D. Practical Integrators 

Finding a closed-form expression for the integral curve 
of some Hamiltonian vector field A is impossible in most 
cases as there is no closed-form solution of Hamilton's 
equations. However, there are some special cases where 
we can find such a solution. 

For example, suppose that in some local patch of phase 
space with coordinates q and p the fundamental 2-form is 
|j a; = A dp, A is an arbitary 0-form, X is an arbitrary 
vector field on phase space. Then 

dA^ dA^ 
oq op 

; _ . 9 d 

and we have 

8 A 8 A 
dA{X) = ^X, + ^X, = u:{A, X) 



d 



d 



= idqAdp)[A,- + A,-,X,-+X,- 



— AqXp 



ApXq. 



Since X is arbitary we can equate coefficients of Xq and 
Xp to obtain 

^_dA d dA d 
dp dq dq dp 

Let c{t) = {qt,Pt) be the integral curve of A with c(0) — 
{qo,Po), which means that for any 0-form / it must satisfy 
the differential equations 



and 



Pt 



dA 



which are Hamilton's equations if A is the Hamiltonian. 

Now, suppose that A(q,p) = T{p) is only a function 
of the momenta, then T — T'{p) d/dq, and Hamilton's 
equations reduce to the pair qt = T'{p) and pt = 
of first-order differential equations with constant coeffi- 
cients, with the solution that the momentum is constant, 
Pt = Po, and qt = qo + T'{po)t grows linearly in t. The 
case where A{q,p) = S{q) is analogous. If we have a 
function A{q,p) — H(q,p) ~ T{p) + S'(g), perhaps the 
Hamiltonian itself, that can be decomposed into the sum 
of a kinetic energy and a potential energy then we can 
easily integrate either term separately, and we can use a 
symplectic integrator to approximate the integral curves 
of H itself. 

In fact we have established a stronger result, namely 
we can find the exact integral curves of a shadow Hamil- 
tonian iJg that differs from H by terms of 0(e) in closed 
form. A symplectic integrator thus not only exactly pre- 
serves the symplectic structure but also conserves the 
value of H (the energy) up to order e for arbitrarily long 
times: unfortunately the integral curves of H and 
usually diverge from each other after a relatively short 
time despite this. This happens even it their equations 
of motion are not chaotic: symplectic integrators are very 
good at conserving energy and phase space volume, but 
they are not particularly good in finding the correct tra- 
jectory through phase space. 

For HMC applications where we only care about ex- 
act reversibility, exact area-preservation, and good en- 
ergy conservation we see that symmetric symplectic in- 
tegrators meet all the requirements, and the divergence 
of the shadow integral curves from the true ones is unim- 
portant. 

Given the fundamental 2-form u} = dq /\ dp we may 
evaluate the Poisson bracket of two arbitrary 0-forms A 
and B, namely 



{A,B} = -u:iA,B) 

'dA d 
dp dq 
dA dB dA dB 
dp dq dq dp 



- {dq A dp) 



dA d dB d 
dq dp' dp dq 



dB_ d_ 

dq dp 



For the Hamiltonian B{q,p) — T{p) + S{q) any integra- 
tor constructed from e^*^ and e^-^ steps will conserve a 
shadow whose BCH expansion may be expressed in terms 
of the Poisson brackets 



(A/)oc=-(/oc), 



We can always find coordinates for whicii tliis is true according 
to Darboux's tfieorem. 



{s,T} = -s"r' 

{S,{S,T}} = -S^ 



,d{S,T} ^ ^,2^„ 
dp 



{T,{S,T}}^T 



= S"T" 
,d{S,T} _ _^„^,2 



dq 
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and so forth. 

For example 



the leapfrog integrator 

exp(i5rS') exp(5rT) exp(i(SrS') is the simplest 

symmetric symplectic integrator (there is a variant in 
which S and T are interchanged). From we find that 
it conserves the shadow Hamiltonian 

E^T^S-^ ({^, {S. T}} + 2{r, {5, T}}) + 0(6t^) 



5t^ 



= H (S"T" - 2S"T") + OiSr"). 

24 



E. Higher Order Integrators 

Let us briefly give some simple examples of more com- 
plicated integrators. The second order minimum norm 
integrator [0-0 is 



t/Sr 



exp {X6tS] exp ( ^^rf ) exp ( (1 - 2\)5tS 



exp {\5tT^ exp [\5tS 



with shadow 



H = T + S + 6t^ 



6A^ - 6A + 1 
12 

1 - 6A 



{S,{S,T}) 

{T,{S,T]]\+0{5t^] 



24 



and it has the free parameter A as well as the integration 
step size St. 

It is interesting to note that if as is usual the ki- 
netic energy is quadratic, T{p) = ip^, then the Pois- 
son bracket {S*, {S*, T}} — S"^ is independent of the mo- 
mentum p, and thus we can find the integral curve of 
its Hamihonian vector field {S, {S,T}} = -2S'S"d/dp. 

The corresponding integrator step e^''^'^'''^'^'-^^^ is called a 
force- gradient integrator step, because it involves second 
derivatives of the potential S. 

We can use the force-gradient step to define a force- 
gradient integrator 



exp(i(5rS) exp(i(5rT) 

X exp 48 StS - 6t^{S, {S,T}} 



: exp(i(5rT) exp(i,5rS') 



with shadow 
H = T + S 



155520 



41{5,{5,{^,{^,r}}}}\ 
+36{{S,T},{S,{S,T}}} 
+72{{S,T},{T,{S,T}}} 
+84{T,{S,{S,{S,T}}}} 
-126{T, {T,{S,{S,T}}}} 

54{r,{T,{r,{5,r}}}}/ 



0(67 



V + 



where we have chosen the integrator parameters to elim- 
inate all terms of 0{5t^) in the shadow. The Poisson 
bracket {5, {S*, {5, T}}} = so the first and fourth Pois- 
son brackets in (||) are also identically zero, however for- 
mula (^) is valid more generally. Note that the middle 
step has combined the Hamiltonian vector fields S and 
{S, {S, T}} because they commute. 

There is no compelling reason to choose the param- 
eters to eliminate the St'^ errors: in general we should 
introduce some parameters constrained only by the con- 
ditions that the leading order term in the shadow should 
be the original Hamiltonian and that the total step size 
should be St, and then adjust these parameters to min- 
imize the cost of our integrator for the specific problem 
it is being applied to. On the other hand, we can build 
integrators whose leading error is St"^ (or (5t^" for any n 
for that matter), without requiring force-gradient steps. 
Nevertheless, integrators with force- gradient steps may 
be cheaper than those without: it would be surprising 
if the optimal coefficient of the force-gradient term was 
exactly zero. 

In HMC for lattice field theory H and H are extensive 
quantities, that is they are proportional to the lattice 
volume V for sufficiently large V, so the leading error is 
proportional to VSt"^"" ii H — H = 0(i5t^"). In order 
to keep the Monte Carlo acceptance rate fixed we there- 
fore need to vary St oc and as the cost VI/St 
of a trajectory of length t is proportional to the num- 
ber of steps and the volume, we may estimate that the 
cost varies as Of course there are many other 
contributions to the cost that have been ignored, but for 
large enough V this suggests that we want to increase n. 



IV. HAMILTONIAN MECHANICS ON LIE 
GROUPS 

A. Fundamental 2-Form on a Lie Group 



The cotangent bundle T*Q over any manifold Q has a 
natural symplectic structure. For the case where ^ is a 
Lie group a point in T*Q may be written as ((?,p) where 
-1 t/Sr 9 & Q p e T*Q[g) is called the momentum or Liou- 
ville form. As explained in Appendix ^ the vectors in 
tangent space at the identity TQ (I) correspond to the Lie 
algebra of left- invariant vector fields on Q, and their 
dual 1-forms 9^ satisfy the Maurer-Cartan equations. 
The momentum may be written in the Maurer-Cartan 
basis as p = Pi6^, where p{ej) = pi6^{ej) = piS^, = pj. 
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We shall choose the fundamental 2-form to be 

ij= dp = -d{p,0'), (6) 

and using the Maurer-Cartan equations it may be writ- 
ten as 

u) = e' Adp, + \p,c)ke^ A6^. 

If F is a 0-form on the cotangent bundle T*Q then 
the corresponding Hamiltonian vector field F — F^ei + 
Fid/dpi in TT*g is defined hy dF = ipu), or dF{y) = 
(^{F, y) for all vector fields y — y^Bi+yid/dpi. Expanding 
this expression gives 

dF 

dF{y) = yF = e,{F)f + —y, 

OPt 

- y) - F% ~ y' F, + p,c]kF^ y\ 

so equating the coefficients of and iji we find (F) = 
-Ft+PjciiF'' and dF/dpi = F\ We thus find that the 
vector field F is 

^ OF ( , dF d 

From this we can evaluate the Poisson bracket of two 
arbitrary Hamiltonian vector fields corresponding to 0- 
forms F and G, 



{F,G} = -u}{F,G) 

dF dG dF 



dG 



C. Poisson Brackets of 5 and T 

We may compute the Poisson brackets of S and T 
from ds) 



{5,T} 

{T,{5,r}} 
{r,{5,{5,T}}} 
{5,{5,{5,r}}} 
{T,{r,{5,T}}} 

{r,{T,{5,{5,r}}}} 



e\S)e,{S) 
-pYeiBjiS) 
2p'e,ej{S)e^{S) 


-pV/e,ejefc(S') 

2pYe^e,ek{S)e\S) 

+2pYe,ek{S)e,e''{S) 



(11) 
(12) 



{{S,T}, {T, {S,T}}} = cyVe'iS) [efce,(5) + e,e,(5)] 



e\S)eke,e,{S) 
-e,e''{S)eke,{S) 



{T,{S,{SAS,T}}}} = 
{{5, T}, {5, {5, T}}} = -2e^(5)eJ (5)e,e, (5) 
{T,{T,{TAS,T}}}} 
{S,{S,{S,{S,T}}}} 



-p'p> p^p^ eiejekei,{S) 



= 



(13) 



Observe that according to equation (|12[ ) {5, {S', T}} d oes 
not depend on the momentum, so just as in §IIIE we 
can use it to define a force-gradient integrator step cor- 
responding to the Hamiltonian vector field 



{e'{S)e,{S) 



d_ 

dpi 



(14) 



Hamiltonian Vector Fields for T and S 



D. Representations 



For HMC we may take the Hamiltonian to be of the 
form H = T + S where the kinetic energy T : r*(y — > R 
is a function only of the momenta which we may choose 
to be of the form 



T=\{p,p) = Wp^ 



(9) 



using the Cartan-Killing metric (§B5). Hence dT/dpi — 
p^, and the potential energy S* : ^ R is a function only 
of the group parameters. 

For the kinetic and potential energy 0-forms the cor- 
responding vector fields are thus 



d 



d 



f = /e, + = p'e, and S = -e,(5)-^ (10) 

dpi dpi 



using (^, where we have made use of the total antisym- 
metry of the structure constants for a semisimple Lie 
algebra, cl^Pjp'' = CjkiP'p^ = 0. 



If [/ : ^ — > Gl(n, (D) = Aut (D is a matrix represen- 
tation of Q then it satisfies U{gh) = U{g)U{h) for all 
g,h € g. We may view any matrix element Uab of the 
representation as a complex valued 0-form as it is well- 
defined over the entire group manifold. The left action 



L 



g . h 1-^ gh induces the map Lg 



Uab 



UabO Lg ac- 



cording to the definition given in §A4, so {LgJJab){h) = 
Uabigh) = [Uig)Uih)]ab = ELi Uacig)Ucb{h) for all 
h, or equivalently LgJJab = Z)c=i Uac{g)Ucb- In other 
words the map Lg^ takes the 0-form Uab to a linear com- 
bination of 0-forms Ucb with coefficients Uac{g) G C. We 
can express this more succinctly by considering U to be 
a matrix-valued 0-form, whence LgJJ = U{g)U. 

Application of the vector field to U gives a matrix- 
valued 0-form BiU whose value at some point g € G 
is eiU{g) — Lg^eiU{T). e.i is left-invariant L*gei — Bi, 



so we have Lg^BiU = Lg,^L*eiU = Lg^Lg- 



e,.LgU 



eiLgJJ 



eiU{g)U = U{g)eiU. This allows us to evalu- 
ate Bill at any point g in terms of the value of eiU at the 
identity. Defining the generators of the representation as 
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Ti = eiC/(I), we obtain e,U{g) = ?7(g)e,[/(I) = U{g)Ti 
or more succinctly eiU = UTi. 

As on the one hand [ei,ej]U = cfjEkU = c'^jUTk, 
and on the other [ei,ej]U = BiBjU — ejeiU = eiUTj — 
BjUT, = UTiTj - UTjT, = U[T„ Tj], we see that the gen- 
erators must satisfy the commutation relations [Ti,Tj] = 
cfjTk upon multiplying on the left by U~^. 

Unfortunately the usual convention ^ is that the 
derivative of a link variable is 

e,U = ~T,U, (15) 

and this is used in most computer implementations. 
This arises from considering right-invariant vector fields. 
Briefly, the right action on a group is defined by Rg : 
h ^ hg, and the induced maps \>y RgJJ — U o Rg 
and R*g^i — Rg-^^^iRgsf- If we assume that e-i is right- 
invariant then it satisfies i?*ei = e^, and following an 
argument completely analogous to that in the text we 
find eiU{g) ~ Rg^eiU{l) since g = Rgl = Lgl and 
Rg^BiU = eiUU{g). We then have to define the genera- 
tors by EiUiJ) — —Ti, leading to eiU ~ —Till. We must 
include the minus sign in the definition of the generators 
for right-invariant vector fields satisfying [e^, ej] — c^jSk 
as otherwise they would not satisfy the commutation re- 
lations [Ti,Tj] = c^jTk. In fact, the usual convention 
erroneously omits the minus sign, but as the commu- 
tation relations are used to derive the Maurer-Cartan 
equations, and thus our fundamental 2-form, the sign is 
significant when computing high order Poisson brackets. 



where we have introduced the quantity Fi = ei{S)T^ 
(q.v., equation ([TtI)). The solution of these equations for 
separate U and P updates (i.e., for a symplectic integra- 
tor) are 

U{t) = cxp{-Pt) U{0) and P{t) = P(0) - tF^. 

The equations of motion for the force-gradient Hamil- 
tonian vector field of ( |l4[ ) is 

P = {S, {S,T}}P = -e,[e,{S)e^{S)) |^ 

= -2e,ej{S)e^{S)T' = -G (16) 

with G = e,ej{S)e'{S)T^ (q.v., equation (|l^), since 
[e„ ej]{S)ei{S) = Ck,je'^{S)ei (S) = 0. 

V. POISSON BRACKETS IN SU{N) 

In order to compute the Poisson brackets it is useful 
to express them in terms of the following set of matrices 
that are in the representation of the Lie algebra 





P 






Fi 




F2 = 


PFi 




F3 = 


V^Fi 


= p''ekp'eMS)T' 


Fa = 


V^Fi 




G = 


J^iFi 


= ei{S)e,e,{S)T' 



(17) 



E. Equations of Motion 

The equations of motion arc most naturally expressed 
in terms of Lie derivatives (§A5). The Lie derivative 
CuT of a tensor field T is its derivative along the integral 
curves of the vector field v, and the definition of the 



Lie derivative given in (A4), (|A5|), and (A6) implicitly 
provides the differential equations defining these integral 
curves. If 1; = ff is the Hamiltonian vector field for 
the Hamiltonian function then these are just Hamilton's 
equations, and we will write T = CfjT. 

For the case of matrix representations we consider the 
matrix elements to be 0-forms as we did in §1VD so we 



may use equation (A4) to obtain U = Cj^U — HU and 
P = CjjP = HP where C/ is a matrix representation 
of an element of Q and P = p^Ti the corresponding ma- 
trix representation of the momentum in the Lie algebra. 
Taking 



H = T + S = p'e,-e,{S) 



d_ 

dpi 

with the explicit forms from (|lo|) , and using the relation 
e,{U) ^ -TiU of (|l|), we find 

U = fU = p'ei{U) = ~p%U = -PU 

P^SP^ -e,{S)^ - -e,{S)r = -Fi 
opt 



Pi 
e^iS) 
p>ejei{S) 
p^e},p>ejei{S) 
p^e£p'^ekP'ejei{S) 



tT{PTi}/a 

tiiFiT,)/a 

lT{F2T,)/a 

tTiF3T,)/a 

ti{FiTi)/a 

tT{GT,)/a 



where V = p^Bi and Fi — e^{S)ei are vector fields (linear 
differential operators) corresponding to the matrices P 
and Fi respectively. For a lattice field theory P, Fi,G, . . . 
will also be lattice vectors, so shall call these quantities 
basic lattice vectors. 

To derive more explicit expressions for the desired Pois- 
son brackets it is useful to use the following identities 
that hold for the fundamental representation of the su(Af) 
Lie algebra ^ for arbitary N x N matrices X,Y, Z, and 



We choose to normalize the traceless anti-Hcrmitian generators 
Ti of the fundamental representation by tr(TiTj) = aSij, where 
a is an arbitrary (negative) constant. For su(3) the Hermitian 
Gell-Mann matrices satisfy tr(AiAj) = 25ij, so our choice 

corresponds to Ti = ^J—apl iXi. Moreover, our definition of the 
kinetic energy is T = ^PiP^ = tr(P^)/2a, and as we observed 
in the introduction changing this normalization corresponds to 
a scaling of molecular dynamics time. One must be careful to 
take all these factors into account when comparing computations 
using different conventions. 
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T, e su(iV) 

4, tr(Xr^) tr(yT'=) = a tr ([X, y]r ) ; (18) 



tr{XT,) tT{YT') = a 



ti{XY) - — trXtrF 



tr[X, y] = tiiXY - YX) = 0; 



(19) 



(20) 



and 



tr i[X, Y]Z) ^ tr{XYZ - YXZ) tr{XYZ - XZY) 
= tr {X[Y,Z]) =ii {[Y,Z]X) 
= tr([Z,X]y); (21) 

from which it follows that 

tr{[X,Y]X)^tY{[X,X]Y) = Q (22) 

and 

Cijk tr(XT*) iv{YT^) iv{ZT^) = tr {[X, Y]Z) . (23) 
Using (^ and (H^) we easily see that 

c%p'p'e^{S)e,ek{S) = ^^(^r'O tr(^^ir^) tr(F2T'=) 



-tr([Fi. 
a V 



(24) 



and as (|l^) leads to 



/[efe,e,](5) = /4£e.(5) = — c^, tr(PrO tr(i^ir ) 



ii {[P,F^]Tk) (25) 



we find using (^ that 

4.pye^(5)[efc,e,](5) 

= ^3 c,,fe tr(PT^) tr(^^iT^) tr ([P, Pi]r^) 

= itr([Fi,P]'^ 



(26) 



Combining equations (E4) and (26) we obtain 



= itr (2[Fi,F2]P+[i^i,P]' 
a V 

We may also deduce from (E3) that 



p*efce,(5) = -tr((P2 - [Fi,P])rfc), 
a 



(27) 



and hence 
P 



p>e''e,{S)ekej{S) = ^ tr ((f2 - [Pi,P] 



(28) 



and 

pV-'e,e^(5)efce,(5) = ^ tr (P| - P2[Pi,P]) • (29) 
From the identity 

we deduce that 

pVe'=(5)efce,e,(^) 

= 4e^-(^V4'y + 24pVe,e,(5)e'=(5) 



-pVe,e,efc(5)e'=(5) 



1 



^4, tr(PiT") iT{PT)ct,m tr(Pr^) tr(Pir'") 
+ ^Cfc,, tr(Pr ) tr(P2T0 tr(Pir'=) 



-^tr(P3rfe)tr(PiT'=) 



-- tr ([Pi, P]' + 2[Pi, P2]P - P1P3 



(30) 



We thus obtain the following expressions for the de- 
sired Poisson brackets 

{5,T} = -tr(PiP)/a 
{5,{5,r}} = tr(Pf)/a 
{r,{5,T}} = -tr(P2P)/a 
{T,{5,{5,r}}} = 2tr(PiP2)/a 
{5,{5,{5,P}}} = 
{T, {T,{5,r}}}--tr(P3P)/a 
{T, {T, {5, {5, T}}}} = 2 {tr(PiP3) + tr(P|)} /a 

{{S, T), {P, {5, P}}} = - tr(3[Pi, P2]P + [Pi, P]' 

-PiP3 + 2P|)/a 

using (1^), (H), (H), and (|^ 
{T,{5,{5,{5,P}}}} = 
{{5, P}, {5, {5, P}}} = -2 tr(PiGi)/a 
{P, {P, {P,{5,P}}}}=-tr(P4P)/a 
{^,{^,{^,{^,P}}}} = 0. 

VI. BASIC LATTICE VECTORS AND TOWERS 

A. Single Link Updates 

We now consider how to evaluate the basic lattice vec- 
tors of (|l7|). This is particularly simple to do in the case 
where there is only a single link variable P, or on a lattice 
if we choose to only update a single link by setting the 
coefficient of the kinetic energy to zero everywhere else as 
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described in § [ A . In this case the potential is of the form ^ 
S = Retr([/X) where X is some constant N x N matrix, 
which in general is neither in the group nor its algebra. 
On a lattice where we are only updating a single link X 
is constructed out of products of other link variables, 
which are themselves constant in molecular dynamics 
time. We find Fi = e^{S)T' = Reti {e,{U)X)T' = 
-RetT{T,UX)T' = -RetT{UXTi)T' = -aTiUX) 
where 7" projects onto the Lie algebra, i.e., the trace- 
less anti-Hermitian part for su(A^). Likewise, F2 = 

VFi = p>ejei{S)T' = RetT(^p> eje,{U)X^T' = 

RetT(p>ej{^T,U)X^T' = -Retrir.pi ej{U)X^T' = 

RetT{T,p>T^UX)T' = Retr(P[/Xr,)T* = aTiPUX), 
and so forth for the remaining quantities in 



Fi = 

VFi = 
F3 = r^Fi - 
F4 = V^Fi = 
G = J-iFi = 



-RetT{UXTi)T' 
RctiiPU XT,)T' 
RetT{P^UXT,)T' 
RetT{P^UXTi)T' 
RetT{FiUXT,)T' 



-aTiUX), 
aTiPUX), 
-aTiP^UX), 
aT{P^UX), 
ar{FiUX). 



B. Lattice Updates 



When we have many links we trivially generalize the 
definition of the fundamental 2-form (|^) to become sums 
over all links 

a; = - ^ dp{i) = -J2 d(pr{i)e\i) 

I f. 



We can compress the notation by letting indices such 
as i also range over all links: that is i — > and the 

implicit sum over the basis of the Lie algebra becomes 
an implicit double sum Of course, we also need 



to augment the structure constants c^j —r ^(^i i i ) = 
c'^^S^/Si'' since the Maurer-Cartan equations do not mix 
links. Similarly, the kinetic energy (]9D becomes 

T = i ^ c{e) {p{£),p{e)) = i c{£)g,,pH£)p^{£) 



^J2c{£)p,i£)p\£) 



where, as discussed in §[A, it is convenient to introduce 
a separate coefficient c{£) in the kinetic energy for each 



We consider the case where the action is Unear in U without loss 
of generality, because if it occurs multiple times we can transform 
it into a form linear in its tensor product, which can be reduced 
into a sum of irreducible representations. For example, the action 
S = RetriUXUX') = Retr[(t/ ® U)X"] where {U r/')ij,fe« = 



UikUjt and X'^^ ^. 



■■ X 



kj 



X'g^ are N'^ X A'"^ matrices, and U ®U 



link. We can extend our compressed notation by implic- 
itly associating a factor of c{£) with each occurence of the 
augmented Cartan-Killing metric, 9{i.ii){j.ij} = 

c{£i)gijSi.ej and hence with every contracted index i. 
With these conventions the definition looks like (||) and 
again. The sums propagate to the Poisson brack- 
ets where the implicit sums over the indices in equa- 
tions (P|)-(p^ also become sums over all links, although 
second derivatives such as 6^6^(5') have bounded support 
for an ultralocal action. It is important to note that the 
implicit factor of q. associated with contracted indices 
means that even though {S*, {S*, T}} does not depend on 
any momentum it still has a factor of c{£) associated with 
each term. If we set c{£') — 6w then only link £ will ap- 
pear in equations (^2|) and (^6|), and the force- gradient 
integrator will therefore only act on that link. 



C. Towers 

The situation would seem to be much more difficult 
when we want to update all of the link variables simulta- 
neously; derivatives like e^^ . . . e^^^ {S) depend on k links 
and it might appear that it will be prohibitively expen- 
sive to compute them. Fortunately we can avoid this 
combinatorial explosion; the key observation is that all 
the Poisson brackets and forces only depend on the basic 
lattice vectors, and these have only a single free lattice 
index. To make use of this we introduce towers of basic 
lattice vectors: a tower T{A, B) is a an array of basic 
lattice vectors T{A, B)i = A^B where yl is a basic lattice 
vector, A is the vector field associated with it, i3 is a sum 
of products of gauge links, and the index i G {0, . . . , n— 1} 
where we call n the height of the tower. 

The basic lattice vectors in (^7|) may be constructed 
from the two towers T(P, B) and T{Fi,B) of heights four 
and two, where B is the stencil of the action S. The 
stencil is the collection of all paths in the action that 
start with a given link. For example, in the case of lattice 
gauge theory without dynamical fermions the action is a 
sum of Wilson loops, each Wilson loop being the trace 
of the product of gauge links around a closed loop. This 
means we can write the action as S — RetT{Ue\~\) + 
Sq where the staple fi is the sum of products of gauge 
links along paths connecting the end of the link £ to its 
beginning, and Sq is independent of Ui, as in §IB. The 
stencil in this case is U^n. This is familiar from the 
computation of the force acting on Ue 



can be reduced into as sum of two irreducible representations 
acting on vectors of dimensions ^N(N — 1) and ^N{N + 1). 



Fi {£) = e,{S)T'^ e, (rc tr [/, n ) 

= Retr(e,(t/^)n)r = Reti{-T,U(n)T 

= -ReiT{UtnT^)T ^ ^ar{Uin). (31) 

The thing to notice here is that we are computing the 
force on the gauge link Ug so the index i is really the 
pair and thus ei{Ut) = for any other link £' ^ £: 

in particular, ei(n) = 0, ei{So) = 0, and ei{Ui\~\) = 
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£i{Ue) n ■ Naturally, we want to compute the force acting 
on every gauge link, and so the stencil computation of 
( ^ ) must be carried out separately for each link: these 
computations can be done in parallel if desired. 

In order to compute the basic lattice vector A' Fi = 
A^ei{S)T'^ we proceed as follows: 

AFi{e) = Ae,{S)T = A RetT{^T,Uin)T' 

= -Retr(T,T(A,{/,n)j)r 
= -aTiTiA,Uin)j 



This is easy to do if we can compute the tower T{A, Uin) 
on the stencil Ui □ • 



D. Algebra of Towers 

It is simple to construct the tower T{A, B) when _B is a 
single gauge link U; we have T{A, U)j = AW = {-A)W. 
This follows from the definitions T{A^ U)o — U and 
A = a'Cj where A = a^Ti, so by induction T{A, U)j+i = 

A+^U = A{AU) = A{-Ayu = a^e,[{-Ayu^ = 

{-Aya'e,{U) ^ {-Aya'i-TiU) = {-Ay+'^U. Indeed, 
this corresponds to a convenient recursive way of con- 
structing the tower, T{A,U)j+i = i-A)T{A,U)j- 

If B is the product ^ of two stencils Bi ■B2 then we may 
use the Leibniz rule for the derivation A, A{Bi ■ B2) = 
ABi ■ B2 + Bi ■ AB2 , or more generally 



AiB^ ■B2)^J2 ( fcj-^'^i ■ -^'"'^^ 



k=0 



The tower on the product Bi ■ B2 is thus the product 
of the tower on Bi with that on B2, T{A,Bi ■ B2) = 
T{A, Bi) ■ T{A, B2), where the product is defined by 



(t(ASi)-T(A,B2)) . = [l)nAB,)k-TiA,B2)j-k 



k=0 



The tower on the sum of two stencils Bi + B2 is even 
simpler, since A{Bi + B2) ^ ABi + AB2- We just have 

T{A, B1+B2) = T{A, Bi) + T{A, B2) where (r(A, Bi) + 
T{A,B2)] =T{A,Bi),+T{A,B2)j. 



E. Pseudofermion Towers 

The principal advantage of updating all links si- 
multaneously is when we include the effects of 
(pseudo)fermions in the dynamics. As described in §[B 
this entails solving a large linear system to obtain the 
quantity X = A4~^(j) needed compute the force {A4 be- 
ing a lattice Dirac operator) and it is worthwhile to reuse 
this solution to update many links. 

We therefore need to compute towers for stencils that 
include outer products such as X (E> X'' . This may be 
done by computing the tower T{A,X) on X ~ J^^^<j). 
Observe that A(j) = as the pseudofermion lattice 
(site) vector (j) does not depend on U — we want to 
follow the molecular dynamics evolution of the gauge 
links and momenta in the presence of a fixed pseudo- 
fermion background. Using the Leibniz rule we get 
= Aict)) = A{MM-^cl)) = AiM)M-^(t)+MA{M-^(t)) 
so AiM-'^(l>) = -M-^AiM)M^^(l). To use this for a 
tower of arbitrary height we generalize this to 



O^A{MM-^(t. 

= MA{M-^4 
for j > 0, and thus 
A{M-^(I)) 



k=0 

E 

k=0 



A-''{M)A''iM- 



k=0 



This translates into the following recursive definition for 
the tower on X 

-1^ 



T{A,X)o = M 
T{A,X)j = -M-^^ 



k=0 



T{A,M)j-kT{A,X)k 



in terms of the tower T{A,Ai) which we already know 
how to compute. Note that we require exactly n inverses 
to construct such a tower of height of height n. 

Yin has suggested an ingeneous way of perform- 
ing a force-gradient update by computing the force twice. 
We should not be surprised that the force-gradient up- 
date e^^ {s,{s,T}} computed out of e^'^^ and e^"^"^ 
steps: recall that according to the BCH formula the com- 



mutator C(e^, e^) 



e[^'^l+ ", hence 



C 



JtS 



C{e 



StS ^StT- 



^-StS^-StT^-StS^StT^StS^-StT^StS^StT 



C 



Jr^[S,[S,T\] + 0{5r'') 



Jr^{S,{S,T}}+0(Sr^) 



^ Here we use the symbol ■ to emphasise multiplication operations. 
Elsewhere we use juxtaposition to indicate multiplication. 
The symbol ■ on the left denotes multiplication of towers, whereas 
on the right it denotes matrix multiplication. 



It is interesting that this can be reduced to only requiring 
two inverses in the case where T is quadratic. There 
does not seem to be a way of using this trick to evaluate 
Poisson brackets, however. 
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VII. CONCLUSIONS 

We have given a formalism for computing integrators 
and the corresponding shadow Hamihonians for lattice 
gauge theories, and we have presented explicit formulae 
for the Poisson brackets up to fourth order and for the 
force-gradient update step. We have shown how to ex- 
press these quantities in terms of basic lattice vectors 
taking their values in the representation of the Lie alge- 
bra, as is needed for the usual formulation of lattice gauge 
theories, and explained how these may be computed us- 
ing towers. The implementation of towers is straightfor- 
ward, as it just requires the substitution of the algebra 
of towers for that of the matrices already used in com- 
puting the force term. The stencils for any action are 
unchanged, and the method is readily applied to pseud- 
ofermions, smeared actions, and so forth. The rules for 
addition, multiplication, and "inversion" of towers are 
given in a recursive form that is easy to implement (al- 
though a recursive implementation is not necessary) . 
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Appendix A: Differential Forms 
1. Differential Forms and Wedge Products 

For convenience we give the definition of a few basic 
operations on differential forms. In some local basis q : 
M ^MJ^ a. fc-form fl e A'' has components 

n= ni^,,j^dq^' A ■ ■ ■ Adq^" 

i<h<-<h<k 

1 ^ 



ii,...,ifc = l 



fc! 

^n,,,„,^dq'^ A--- Adq'" 

^ ^ n^.^.^^dq""' A--- Adq^' 

{n^,...^,dq-^ A---Adq--)^^S^ 



where Sk is the symmetric group acting on 1, . . . , fc, and 
(• • ■)s^ indicates the average over elements of the sym- 
metric group. The wedge product satisfies 

aA/3 = {-lY^'l3AoL ol e h-'' , (3 e k^' 

Antisymmetry; 
aA/3A7 = aA(/3A7) = (Q:A/3)A7 

Associativity. 

In terms of the components in local coordinates this 
means that |^ 



a.Aj3 = (a^^...nJn^^,...n,:,,dq'''A ■ • • A dq^''-^"' 



k + k' 



(fc + fc') 



^ ^ a^,....J.,^,....,^„dq-^A---Adq^>'+^'. 



2. Exterior Derivatives 

The exterior derivative d : ^ A'^'^^ is a linear an- 
tiderivation, so 

d{a + f3) = da + df3 Linearity; 
d{aA(3) = (da) A 13 +{-!)'' a A dfB a e A'' 

Anti-Leibniz; 

d^a = 

dF{x) ^ xF F e AO. 

The exterior derivative dF for a 0-form F is defined to 
be dF{x) = xF for any vector field x: if we evaluate this 
in a local coordinate system we find that 

dF{x) = xF= (x^^) F = (^) x^ 



so 



dF 

dF = ^dq\ 



Likewise, in a local coordinate system the exterior deriva- 
tive of a fc-form $7 e A*^ is 

dn = d A • • • A dq'" 

■^dq^ Adq'' A--- Adq'". 



fc! dqi 



This follows from the anti-Leibniz rule d{a(3) = da A 
(3 + ad(3 applied to the case where a = ^ii...ik G f^^d 
(3 = dq^^ A - ■ ■Adq^'' because the second term vanishes (by 
induction on fc) using the condition = for the basis 



Our convention is that each independent component occurs once 
in the sum: another convention is that each such component 
occurs fc! times — once for each permutation of its indices. 



For the other convention the numerical coefficient in this formula 
is l/(fc!fc'!): caveat emptor. 
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forms which are exterior derivatives of the coordinates 
q\ <fq' = . 

In particular, for a 1-form £ we have 



dB 



09, 



dq^ A dq\ 



so applying the 2-forni dO to two arbitrary vector fields 
X and y gives 

do 

d9{x,y) = g^i^'y' - x'y^) 



d 



d 



- dqJ ^^'^^ ~'' 'dqi " dqi ' ' " " dq^ 

= xe{y)-yeix)-0, [x{f)-yix')] 
= xe{y)-ye{x)-e{[x,y]). (Al) 

This provides an elegant coordinate-independent defini- 
tion of dO in terms of the commutator of the vector fields 



d , d 



d , d 



[x,y]^xy-yx^x^—y^—-y^—x^— (A2) 



dq^ dqi 



dq- ^dq-JdqJ^^^ ^ ' dq^ dqJ ' 

which is itself a vector field since the last term involving 
second derivatives vanishes by symmetry. Note that if 6 
is exact, that is = dF, then the identity d^F{x,y) ^ 
xdF{y) — ydF{x) — dF{[x, y]) — xyF — yxF— [x, y]F = 
holds automatically. 

For an arbitrary (fc — l)-form G A'''^^ we may derive 
the corresponding identity, 

k 

dn{xi, . . . ,xk) = ^{-iy^'^Xin{xi, . . . ,Xi, . . . ,xk) 
1=1 

— ^ ^ ( 1) '^•^^ f2( [ajj, ajj] , a?! , . . . , , . . . , Xj , . . . , a;^,) 

l<i<j<k 

where x indicate that the variable x is omitted. We 
observe that for k — 3 the invariant expression for the 
exterior derivative is 

du){x, y, z) = xu){y, z) — yuj{x, z) + zijj{x, y) 

-Lc{[x,y\, z) + Lo{[x, z\,y) - uj{[y,z\,x). (A3) 

3. Interior Products 



4. Induced Maps 

If (7 : Al — !■ A^' is a diffeomorphism, then there is a 
natural induced map a* : k°{M') A°{M) defined by 
cr*f -p^ f{<jp) for all / e A"(X') and p e M. This 
map may also be written as cr*/ = / o c, and is called a 
pull-back. Another way of saying this is that the following 
diagram commutes 

M ^ M' 

E, 

If (7^^ exists then there is a corresponding pull-back 
map (cr^^),, and it satisfies the relation ((7^^),cr»/ — 
(ct^^),(/ o a) — f o a o a^^ — /, and thus we see that 
(ct^^), — (ct,)^^, and we may denote both of these un- 
ambiguously as cr^^. 

If a; G TM is a vector field on M then there may 
be a push-through map a* : TM — s> TA^' defined by 
a*x — o a; o CT* if this exists. For any / 6 hP{M') 

x{a,f)\^. The 



and p £ Al this means that a*x{f) 
corresponding commutative diagram is 



-1 

A°{M) 



A°{M') 
A°{M'). 



The existence of the diffeomorphism cr^^iA^'— i>A1isa 
sufficient but not necessary condition for a^^ and hence 
a* to be well-defined. 

We may define further induced maps ^ such as the 
pull-back of one-form fields cr, : A^(Al') A^(Al) as 
cr*0 — o o a* 



TM 



TM' 



AO(A() ^ A"(A('), 

and so forth. 

In the special case where ct : A^ — > A^ is an autodif- 
feomorphism then the push-through maps always exist. 



The interior product i : TM x A''' A'^^^ is the op- 
eration that inserts a vector as the first argument of a 
/c-form to yield a fc — 1-form. It is formally defined by 
the axioms 

ix{cy- + /3) = ix<^ + a, /3 e A*^ Linearity 

i^(aA/3) = i^{a) Af3+i-l)''aAi^p 

a. e A'^ , f3 G A''' Anti-Leibniz 
i^F ^ F eA° 

ia:^{xi, . . .,Xk-l) = fl{x,Xi, . . .,Xk-l) ^ £ A*-' 

so we see that it too is a linear antiderivation. 



5. Lie Derivatives 

Suppose now that we have a smooth one-parameter 
family of diffeomorphisms ct : IR x A^ ^ A^, which we 



One must be careful with the notation introduced here, as there 
are a whole family of mappings that we have given the same 
name, cr, : A'^(AI') — > A''{M) Vfc, and the equation (7,0 = 
at o0 OCT* involves two of them. If we were to call these induced 
mappings on forms crj : A''{Ai') — > A*(A^) then the equation is 
less ambiguous, a^O = a^l o o a* . 
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will also write as at : M. ^ M- Using this map we can 
define a derivative with respect to the parameter which 
is called a Lie derivative. For any 0-form F we define 



dt 



d{F o at) 



dt 



vF 



(A4) 



where v is the linear differential operator — the vector 
field — that is tangential to the curves cr(t,p) passing 
through ct(0,p) = p G >1 at t = 0. 

The Lie derivative of a vector field y G TAA can be 
deduced from the requirement that Cx be a derivation 

Cx{A ®B) = {C^A) ®B + A® C^B 

for any tensors A and B, and that it commutes with 
contractions 

Cx{yF) = {Cxy)F + y{CxF), 

Cx{e{y)) = {CMy) + o{Cxy) 

and so forth. Applying these rules to the 0-form yF 
obtained by applying the vector field y £ TA^ to F e 
k^{M.) we have C^iyF) = xyF and also Cx{yF) = 
{Cxy)F + y{CxF), hence 



{Cxy)F = xyF - yxF 
and, as this holds for all F, 

-Cx?/ = [x,y]. 



(A5) 



We may apply a similar argument to evaluate the 
Lie derivative of a 1-form 6 G A^(A^). On the one 
hand Cx{0{y)) — x9{y), while on the other Cx{0{y)) ~ 
iCxe){y) + e{Cxy), so using (P) 



{Cxe){y) = xe{y) - e{[x, y]) = de{x, y) + ye{x) 
- {ixde){y) + d{e{x)){y) = {ixde){y) + idix9)iy) 
= {ixd + dix)9{y), 

hence 

CxO = {ixd + dix)6. 

This suggests that the Lie derivative of any /c-form may 
be expressed as 

Cx = ixd + dix, (A6) 

and this is indeed the case as the operator ixd + dix is a 
derivation 

{ixd + dix){cx. A/3) 

= ix [{da) A 13 + (-!)'-'« A d0\ 

+d [{ixcx) A (3 + {-l)''a A ixl3] 

= {ixda) A 13 + (-l)'=+i(da) A ixl3 

+{^l)''{ixa) Ad(3 + {-l)^''a A ixdf3 
+{dixa) A/3 + {-l)''-\ixa:) A df3 
+{-l)''{dcx) A ix(3 + {-l)'^^OL A dixl3 

= [{ixd + dix)a] A (3 + a A {ixd + dix)f3 



for all a G A*^ and /3 G A*^ , and for 0- and 1-forms F 
and 6 

CxF ^xF = dF{x) = ixdF + dixF, 
CxO = {ixd + dix)6. 

The second term in the first equation is zero because 
ixF = by definition. 



Appendix B: Lie Groups 
1. Left- Invariant Forms 

A Lie group is a manifold that has a group structure 
defined by C°° multiplication {g, h) M- gh and inverse 
g ^ g^^ operations that satisfy the group axioms 

9{9y) = {99')/' = 99' 9" y9,9',9"^G Associative 
g~^g = gg^^ — I ^9 & G Inverse 

with I being the identity element of the group. If we con- 
sider a point g (z Q as being "fixed" then left multiplica- 
tion by g is an autodiffeomorphism of tj, Lg : g' ^ gg', 
with Lgh = LgoLh by associativity, LgoL^g' — g{hg') = 
{9^)9' — Lghg' for aU g' G Q. Clearly Lg-i — {Lg)~^ too. 

As for any such diffeomorphisms we can define the cor- 
responding pull-back maps on forms and vectors, Lg^F = 



FoLg, L*gV 



Lg-i^ovoLg^, and Lg^9 



L„ oOoL* We 



may use these maps to define left-invariant vector fields 
and forms; for example, a left-invariant 1-form satisfies 
the condition 9 — Lg^B. 



2. Lie Algebra 



The only left-invariant 0-forms are constants, as if 
F = Lg^F (V5 G g) then F{g) = F{Lgl) = Lg^F{I) = 

If u = L*gU and v ~ L*gV are left-invariant vector fields 
in the tangent bundle TQ then their commutator is also a 
vector field, and furthermore it is also left-invariant since 
f^[u,v] = [L*gU,L*gV] = [Lg-i^ouoLg^,Lg-i^ovoLgJ = 
Lg-i^ o [u, v] o Lg^ ~ L*g[u, v]. If a left-invariant vector 
field V vanishes at the identity, v{F)\i = (VF G A°^), 
then it must be identically zero everywhere, as v{F)\g = 
v{F)oLg]^ ^ [LgMF)], - [Lg_^oL*gV{F)]^ = 



v{Lg^F) =[v{FoLg)L = 



Consider a set of left-invariant vector fields {e^} in 
TQ called generators whose values at the origin are lin- 
early independent; any linear combination of the gen- 
erators with left-invariant (constant) coefficients is also 
left-invariant. Conversely any left- invariant vector field u 
must be a linear combination of this type, since its value 



Note that L„-i 
9 » 
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at the origin is u\i = ^i^^ ^ hence 

n — M*ei = everywhere. Left-invariant vector fields 
therefore form a hnear space with constant coefficients. 
In particular, the commutator of any left-invariant vec- 
tor fields must be a linear combination of the generators, 
[si, Ej] — c^jEk where the c*j £ IR are called the structure 
constants. This makes the linear space of left-invariant 
vector fields into a Lie algebra. 

Any left-invariant vector field v has an integral curve ^ 
c : R — ^ satisfying c(0) = I. Along this curve we have 
an Abelian subgroup of Q satisfying c(s + t)— c(s)c(i), 
so it is naturally to call c an exponential map, and write 
it as c{t) = exp{vt). If we view this as a function of v 
then this defines a local flow of v, and is a map from the 
Lie algebra into the Lie group, exp : TQ Q. 

The commutator of two elements g,h G G is defined 
to be C{g.h) = g^^h~^gh; in a neighbourhood of the 
identity where g = cxp(ut) and h = exp(t)i) we have 

C{g, h) — exp{—ut) exp{—vt) exp{ut) exp{vt) 
= (l-ut+ ^{utfj (l-vt + ^{vt)'^ 



X [I + ut+ ^{utf 



I + vt+^{vtf 



Oit^) 



I + [u, v]t^ + 0{t^) = exp([i6, v]t^) + Oit^). 



3. Maurer—Cartan Equations 

The commutation relations may be succinctly ex- 
pressed in terms of the cotangent space T*G. We intro- 
duce a set of left-invariant 1-forms 0* (called a frame or 
repere mobile) dual to the generators 6^{ej) — 5^ From 
(Al) we have 



so expanding the 2-form dO^ = al^„9™ A 9" in terms of 
the basis 2-forms 6™ A 9" we have 

d0(e„efe) = <„0"A0"(e„efe) 

= <„ {0"(e,)0"(e,) - 0™(efc)0"(e,)} 

thus the left-invariant forms 0' satisfy the Maurer- 
Cartan equations d9^ = ^h'^^jk^'' ^ everywhere. 



4. Adjoint Representation 

For any Lie algebra the adjoint representation is de- 
fined by ad(a;)y = [x, y\. This is a representation of the 



Lie algebra because for any z 

[ad(a;), ad(y)]2; = ad(a;) ad(y)2; — ad(y) ad(a;)2; 

= ad(a;)[y,2;] - ad(y)[a;,2;] 

= [a?, [V.z\] - [y, [x,z\] = [[x,y\,z\ 

where we used the Jacobi identity in the penultimate 
step, and thus [ad(a;), ad(y)] ~ ad([a;,y]). In terms of 
basis vectors we have ad(ei)ej = [ei,ej] = c^jSk, giving 
the explicit matrices ad(ei)^ = . 



5. Cartan— Killing Metric 

We may use the adjoint representation to define 
the Cartan-Killing metric on the Lie algebra as a 
trace, {x,y) = tr[ad(a;) ad(y)]/CA where Ca is a con- 
stant; in terms of the basis vectors gij = {ei,ej) = 
tr[ad(ei) ad(ej)]/Cyi = c^g^c^-j,/CA- For a semi-simple 
Lie algebra the Cartan-Killing metric is non-singular and 
has an inverse satisfying g^-'gjk = ^fe- For a simple Lie 
algebra the adjoint representation is irreducible, so by 
Schur's lemma the invariant Cartan-Killing metric is a 
multiple of the unit matrix; we shall choose the con- 
stant Ca such that this multiple is unity. For su(A^) 
where the generators in the defining N dimensional fun- 
damental representation Ti satisfy the commutation re- 
lations [T^,Tj] = c^jTk and are normalized such that 
tvTiTj = aSij the Cartan-Killing metric is explicitly 
gij = Sij with Ca = 2aN. 

For semi-simple Lie algebras we can use the Cartan- 
Killing metric and its inverse to lower and raise indices 
at will, for example we shall define = g^^Pi, and cor- 
respondingly we have an invariant quadratic form for 1- 
forms, {a, (3) = g^^aifHj where a = atO"^ and (3 = PiO"^ . 
We also note that the quantity Cijk = guc^jk — ^Cifcj is to- 
tally antisymmetric, because ([e^, ej], e^) — cjj (e^, e^) = 
cfjgtk = Ckij, and 

Ca{[X,Y],Z) =ti(^&d{[X,Y])&d{Z)' 
= tr([ad(X),ad(y)]ad(Z)^ 
= tr(ad(X) ad(r) ad(Z) - ad(y) ad(X) ad(Z)) 
= tr(ad(Z) ad(X) ad(r) - ad(X) ad(Z) ad(r)) 
= tr([ad(Z),ad(X)]ad(l^)^ 
= tr(ad([Z, X]) ad(r)) = Ca {[Z, X],Y) , 

hence Cijk — Cjki — ^kij- 



strictly speaking this is only true locally: to be precise we should 



write c : X — >■ M where X C R, is a neighbourhood of zero. 
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